function hts = calcAhInv(invP,etaModel,q)

k = 1/invP.disT(1,2);


%Mass matricies
T = sdiag(invP.mesh.Af*(1 + k*invP.tauModel.*(1 - etaModel)));
Tt = sdiag(invP.mesh.Af*(k*invP.tauModel.*(1 - etaModel)));
M = k*invP.mu*speye(invP.mesh.ne);
S = sdiag(invP.mesh.Af*(invP.sigmaModel.*(1 + k*invP.tauModel)));
St = sdiag(invP.mesh.Af*(k*invP.tauModel.*invP.sigmaModel));

A = invP.mesh.CURL'*(S\T)*invP.mesh.CURL + M;

hts = cell(invP.nt,1);
ets = cell(invP.nt,1);
jts = cell(invP.nt,1);

rhs = -invP.mesh.CURL'*(S\q{1});

hts{1} = A\rhs;
ets{1} = (S\T)*invP.mesh.CURL*hts{1} + S\q{1};
jts{1} = (T\S)*ets{1} - T\q{1};

for ii = 2:invP.nt
    rhs = M*hts{ii-1} - invP.mesh.CURL'*(S\St)*ets{ii-1} ...
                        + invP.mesh.CURL'*(S\Tt)*jts{ii-1} ...
                        - invP.mesh.CURL'*(S\q{ii});
    hts{ii} = A\rhs;
    ets{ii} = (S\T)*invP.mesh.CURL*hts{ii} + (S\St)*ets{ii-1} ...
        - (S\Tt)*jts{ii-1}  + S\q{ii};
    jts{ii} = (T\S)*ets{ii} - (T\St)*ets{ii-1} + (T\Tt)*jts{ii-1} - T\q{ii};
    
end
